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Abstract 

The single graphene layer is a novel material consisting of a fiat monolayer of carbon atoms 
packed in a two-dimensional honeycomb-lattice, in which the electron dynamics is governed by the 
Dirac equation. A pseudo-spin phase-space approach based on the Wigner- Wcyl formalism is used 
to describe the transport of electrons in graphene including quantum effects. Our full-quantum 
mechanical representation of the particles reveals itself to be particularly close to the classical 
description of the particle motion. We analyze the Klein tunneling and the correction to the 
total current in graphene induced by this phenomenon. The equations of motion are analytically 
investigated and some numerical tests arc presented. The temporal evolution of the electron-hole 
pairs in the presence of an external electric field and a rigid potential step is investigated. The 
connection of our formalism with the Barry-phase approach is also discussed. 



I. INTRODUCTION 



Graphene can be considered as one of the most splendid functional materials. This has 
been proved by the quick science respond to the novel fascinating experiments performed by 
A. K. Geim and K. S. Novoselov pQ. Graphene represents a single layer of sp^-bonded carbon 
atoms, which are densely packed in form of a benzene ring structure. This ideal planar 
structure has been used to describe properties of many carbon-based materials including 
graphite (that can be viewed as a large number of superposed graphene sheets). This new, 
strictly two-dimensional material displays unusual electronic properties arising from the bi- 
conically shaped form of the Fermi surfaces near the Brillouin zone corners (Dirac points). 
In a quite wide range of energy, electrons and holes propagate as massless Fermions and 
their behavior reproduces the physics of quantum electrodynamics but at the much smaller 
energy scale of the solid state physics. 

For example, a typical superconductivity phenomenon, like the Josephson effect, have 
analogs in the p-n junction. In [2] the particle transport trough an interface between a normal 
and a supercoductor material (N-S) is compared with the analogous inter-band tunneling in 
a p-n graphene junction. It is shown that for excitation energies, which are small compared 
to the superconducting gap, the Dirac Hamiltonian of an p-n junction displays the same 
excitation spectrum as an N-S junction. 

New experiments preformed in graphene-based materials showed the evidence of the 
simultaneous occurrence of relativistic-like and superconducting transport that opens the 
possibility to study the implication of the relativistic-like electron behavior in the solid state. 
In [3] Heersche et al. studied the supercurrent flowing through a simple device constituted 
by two superconducting electrodes on the top of a carbon monolayer. They observed the 
occurrence of some new properties displayed by massless particles in graphene as for example 
the integer quantum Hall effect and the Aharonov-Bohm effect [U E]. Moreover, a large 
Rashba splitting (corresponding to an energy shift of 225 meV) of the tt states in an epitaxial 
graphene layer on a Ni substrate has been reported in p]. 

In particular, evidence has mounted that the scattering of electrons near the Dirac point 
in graphene-superconductor junctions differs from the analogous Andreev scattering process 
in normal metals and that the quasi-ballistic transport in graphene sheets is only weakly af- 
fected by external sources of disorder (defects or impurities). In fact, Dirac fermions are quite 
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immune to the localization effects and it has been observed experimentally that electrons can 
propagate without scattering over distances of the order of micrometers [7j. Furthermore, 
in high-density low-temperature regimes, the mobility is roughly density- and temperature- 
independent. Although the interactions with the underlying substrate are largely responsible 
for the relaxation properties of particles in graphene (possible sources of scattering include 
adsorbents and defects in the graphene lattice, ionized impurities in the silicon oxide sub- 
strate, surface charge traps, interfacial phonons and substrate ripples [8]), the exact nature 
of the scattering that limits the mobility of electrons in graphene devices remains unclear. 

Graphene is a semiconductor, whose band gap is exactly zero and the velocity of the 
charge carriers is over a wide range of energy independent of the momentum. Graphene is 
expected to be in a low-conductivity state when the Fermi energy approaches the Dirac point 
where the density of states vanishes. A gate voltage can, however, modulate the density of 
states in graphene and switch between the low-conductivity state at the Dirac point and 
the high-conductivity states elsewhere. The charge mobility in graphene layers attains large 
values that cannot be reached in conventional semiconductors (mobility of the order of 
lO^cm^V^^s^^ have been recently measured [9j). Because of this high electronic mobility 
and the capability of being tuned from p-type to n-type doping by the application of a gate 
voltage, graphene is an interesting candidate towards possible applications in carbon-based 
electronics devices. In particular, some applications are already devised by various groups, 
for example in designing electronic building blocks [TOl [TT] or spin injection devices [12]. 

Moreover, when the Fermi level approaches the Dirac point, the density of states vanishes 
and it is expected that also the conductivity becomes strictly zero. On the contrary, the 
theoretical prediction of Fradkin, given in [13], concerning the presence of a residual minimal 
charge conductivity, was confirmed by experiments. The main reason of this phenomenon 
concerns the difficulty of localizing Dirac-like particles in a single band. The possibility 
to perform easily band-to-band transitions, provided by the gap-less Dirac-like form of the 
Hamiltonian, reveals that the particles can travel over long distances (or penetrate a poten- 
tial barrier) without creating a reflected component by converting itself in a electron-hole 
excitation. Because of the strong similarity with relativistic quantum mechanics, the tun- 
neling of an electron through an p-n graphene junction, where conduction-like states are 
converted into hole-like states (and viceversa) is denoted as Klein tunneling. It represents 
the tunneling of a particle into the Dirac sea of antiparticles (represented by the almost filled 
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hole band). In several recent experiments, this unusual coupling of electron- and hole-like 
dynamics have been investigated [E]. Klein tunneling gives rise to some unusual behavior 
of the charge transport when the Fermi level approaches the Dirac point, where the valence 
and conduction bands meet. In particular, an unusual interesting transport phenomenon 
of relativistic-like particles concerns the normal incidence of a particle-antiparticle beam on 
a square potential barrier. When the incidence angle becomes equal to zero, the barrier 
becomes completely transparent (Klein paradox). This result is characteristic for the Dirac- 
like dispersion relation of the Hamiltonian and contrasts the electron transport in normal 
(nonrelativistic) devices, where the band-to-band transmission probability is always smaller 
than one. 

In solid state physics, we are typically interested in macroscopic phenomena, which are 
slowly varying in time and smooth in space apart from variations on the atomic scales. The 
language used to describe electron transport is derived from the semi-classical picture of 
the dynamics where the electrons respond to external fields like point particles. There have 
been overwhelming evidences that such a simple picture cannot give complete account of 
first-order effects in the fields. 

The development of efficient quantum computational methods is thus a crucial aspect 
in the study of new devices where quantum-mechanical effects play a dominant role. Dif- 
ferent approaches based on the density matrix, non-equilibrium Green's functions, and the 
Wigner function have been proposed to achieve a full quantum mechanical description of 
the electron transport [15j. Among them, the Wigner-function formalism is the one that 
bears the closest similarities to the classical Boltzmann equation, so that this formalism can 
be considered as a natural choice to derive quantum corrections to the classical phase-space 
motion. Furthermore, a phase-space approach may appear more intuitive compared with the 
more abstract density matrix and Green's function formalism. The phase-space formulation 
of quantum mechanics offers a framework in which quantum phenomena can be described 
with a classical language and the question of the quantum-classical correspondence can be 
directly investigated [T6] . 

For these reasons, an approach where both the kinetic characteristic of the particles and 
the pseudospin degree of freedom are described in a full-quantum mechanical framework, 
seems to be a promising approach to shed light on these particularities of graphene. The 
close similarity between the classical mechanics and a quantum kinetic framework, which 
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characterizes the Wigner single-band formahsm, is generally lost when we address the many- 
band dynamics. In fact, a straightforward extension of the standard definition of the Wigner 
function leads to very complicated multi-band systems, where a one-to-one relationship 
between band and distribution function cannot be found. In general, it is not possible to 
define a quasi-distribution function associated to a single kind of particles (holes or electrons) 
and whose marginal distributions (for example the integral with respect the momentum) 
represents some expectation values of such particles. 

In this contribution, we therefore present a Wigner-like multiband formalism and study 
the effect of Klein tunneling. In sec. |TT] the derivation of the quasi-diagonal equations of 
motion is presented and our approach is compared with some preexisting methods. In sec. 



Ill we discuss the particles motion in the presence of a uniform electric field. Particular 
emphasis is given on the close similarity between the description of the Klein tunneling 
provided by our method and the classical particle transport. The numerical difficulties for 



a direct solution of the transport equations are discussed in sec. IV and an asymptotic 
approach is proposed. Finally, in sec. |V] we study the particle motion in the presence of a 
rigid barrier. 



II. WIGNER FORMALISM FOR THE QUANTUM TRANSPORT IN GRAPHENE 

The atomic structure of graphene is characterized by two types of bonds and exhibits the 
so-called planar sp^ hybridization. The a bonds are strong covalent bonds responsible for 
most of the binding energy and for the elastic properties of the graphene sheet. However, 
since the upper (lower) bound of the a (a*) band is quite faraway the Fermi energy (more 
than 4 eV and 8 eV at the F point for the a and the a* orbital, respectively), bonding and 
anti-bonding a bands can be safely neglected when addressing the electronic properties of 
graphene. The half-filled vr bands are responsible for the charge transport. The first who 
studied the graphene band structure was P. R. Wallace in 1946 by using a tight binding 
approach [T7]. Subsequently, more refined models were derived, providing a reliable theo- 
retical basis for the description of the electronic properties of this material (an exhaustive 
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bibliography concerning this models can be found in [18]). The Hamiltonian [T^ [19] 



reproduces the spectrum of an electron-hole pair in a graphene sheet lying in the x — y 
plane, in the presence of an external potential U{r). Here, vp is the Fermi velocity, cr = 
{ax, Uy^az) denotes the Pauli matrices vector and ctq the identity 2x2 matrix. The valence 
and conduction bands are usually denoted as pseudo-spin components of the particle. 

From a technological point of view, the direct integration of a graphene sheet in a de- 
vice could cause some disadvantages. They are mainly related to the absence of an energy 
gap between the particles and holes at the Fermi surfaces that prevents the electrons to be 
electrostatically confined in graphene. If compared with other open gap carbon-based struc- 
tures, as for example carbon nano-ribbons (where band gaps of nearly 1 eV are observed), 
the absence of a gap allows high current flows also in the off state. This strongly limits the 
application of a graphene sheet as a suitable channel in a carbon-based FET. For this reason, 
we derive our evolution model for a more general Hamiltonian than Eq ([2]), containing an 
energy gap A at p = 0: 



We establish the particle equation of motion in the quantum kinetic formalism by defining 
a suitable multi-component Wigner function. From a technical point of view, one of the main 
purposes of our approach is to describe the particle evolution by a set of Wigner functions 
/ij(r, p) in such way that each function fij is the Wigner transform of a mixture of electronic 
states belonging only to the i-th and j-th band. This ambitious goal would require the 
diagonalization of the pseudo-spinorial Hamiltonian "H of Eq. ([s]) in the momentum as well 
as in the position space. However, this is in general impossible, due to the non-commutativity 
of these operators. In the following, we propose a procedure that tries to define a set of 
basis states that diagonalize "as much as possible" "H. For that reason, the corresponding 
set of Wigner functions will be denoted as a "quasi-diagonal" representation. 

We study the electron-hole pair system by means of the Weyl quantization procedure. 
For the sake of completeness, we recall briefly the basic mathematical tools used in the Weyl 




(1) 




HA=no + Aao. 




(3) 
(4) 
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formulation. Given a differential operator A (defined on a suitable Hilbert space 
function h, the Weyl map W [A] (h) = Ah, is defined as 



and a 



Ah 



It establishes a unique correspondence between A and a function ^(r, p) which is denoted 
as the symbol of the operator. Here d is the dimension of the position and momentum 
space. In the framework of the Weyl quantization procedure, a mixed state is defined by 
the density operator 

S[h] = J p(x,x')/i(x') dx' 



whose kernel is the density matrix. The Weyl symbol 5 = W 



-1 



S 



is obtained by applying 



the inverse of the Weyl transformation (Wigner transformation) to the function p{x, x') [20] : 



p r 



2'" 



e 1 



p-TJ 



(5) 



The von Neumann equation 



ih 



dt 



(6) 



gives the evolution of the density operator and expresses the evolution of the system in op- 
erational form. By using the Weyl operator, this equation can be mapped into an evolution 
equation defined in the phase plane r — p. The symbol associated to the graphene Hamilto- 
nian given in Eq. (|3]) is H (r, p) = W"^ V. = 'Ha(p) + o"of^(r) where Ha = vpcr ■ p + A (Tq 
(in this simple case, the usual quantization —ihV p holds true). We consider the density 
operator S' = Q S Q"^ where B (r, Vr) is a unitary 2x2 matrix operator and the super- 
script t denotes transposition and conjugation. A convenient quantum kinetic description 
of the electron-hole pair motion can be obtained if we exploit the link of B with the symbol 
B(r,p) = W^^ B . In particular, we require that B (r, p) diagonalizes the Hamiltonian 
"H (r, p) locally in the position and in the momentum space. We have 




VE + A e-'^^VE-A 
e'^^VE-A -y/E + A 



(7) 
(8) 
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where A(p) = (T^i?(p), the relativistic-hke spectrum of the graphene sheet E = A/f|,|pp + A^, 



and e*^p = -^y==. Equation (|6|) transforms to 



ih 



dS' 

~dt 



H'S' 



(9) 



where Ti' = Q Ti &^ . By applying the operator W ^ to Eq. (9), we obtain the final equation 
of motion for the symbol S' = VV~^ S' in the phase-space (r, p) (details of calculations are 



given in Appendix VII A): 



^h— = [U' + Aip),S'] 



(10) 



where the brackets denote commutation -B]^ = A-kB — B-kA. The star-Moyal product 
•k is defined as 



(11) 



where the arrows indicate on which operator the gradients act. The symbol W (r, p) is given 

by 



W (r,p) = 0^f/(r)^e"^ 



(12) 



and writes explicitly as 



W'(r, p) = -1-^ 1 (^p + et (^p - t/(r')e^('-^') '^ d^ dr' . 



Equation ( 10 ) is given in terms of the Moyal commutator and defines implicitly a non-local 



evolution operator for the matrix- Wigner function S'. It requires the evaluation of infinite- 
order derivatives with respect to the variables r and p. The commutators appearing in Eq. 



(10) can be written in integral form as 
1 



[A, 5'] 



1 



A(p+^m) 5'(r',p)-5'(r',p)A(p-^M 



(13) 



W ( r - ^77, p + ^/.] S' (r', pO - S' (r', p') W (r + ^^7, p - 



X e^(''^'-') ''+^(P-P')-^ dn dr' d-q dp' . 

(14) 



The commutator of Eq. (13) describes the free motion of the electron-hole pairs in the 
upper and lower conically shaped energy surfaces S^. One of the principal aims of the 
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diagonalization procedure of Eq. ([T]) was to derive an equation of motion, where the free 
motion is described in terms of the evolution of two non-interacting particle populations. 
This is achieved since A is a diagonal matrix. The free evolution of the particles if ) 
belonging to the upper (lower) part of the spectrum is described by 



dt (27r) 



2' 



where we defined the components of the matrix S' as 

5'^(2.fi)-^f{!<''P> ^'''-P' I . (16) 
/'(r,p) /-(r,p) 

These equations describe the free quantum mechanical motion in the band structure defined 



semi-classically by the function -E'(k) = y Vph? + A"^ and generalizes the mass term 
present in the parabolic band approximation. It should be noted that our procedure is 
derived in a full quantum mechanical context, without invoking the usual generalization of 
the semi-classical motion to the quantum mechanical one, where the substitution k — t- — iVr 
in the semiclassical expression of the energy spectrum E{k.) is assumed. As expected from 
a physical point of view, the coupling between the bands arises from the presence of an 
external field U{r) which perturbs the periodic crystal potential. This is described by Eq. 



(14). 



In order to illustrate the main features of the pseudo-potential W(r,p), in fig. [T] we 
depict the explicit form of W when the external potential U (r) (represented in the sub-plot 
[T]-a) is a barrier. Form Eq. (12) we recognize that W is a 2 x 2 matrix depending both 



on the position r and the momentum p. We note that some p-dependent corrections to 
the potential arise around p^; = 0, whereas the pseudo-potential stays practically identical 
to U for greater values of the momentum. We remark that this characteristic reflects the 
presence of a singular behavior of the particle-hole motion in the proximity of the Dirac 
point. It will be addressed in more details in the following sections. In order to highlight 
the modification of the pseudo-potential when the parallel momentum py changes, in fig. 
|2| we depict the component [W(r, p)]_|_+ (in a single band description of the dynamics, it 
represents the potential "seen" by the particles in the S+ band) for different Py. The plot 
shows that for large values of Py, the in-band component of the pseudo-potential coincides 
with the external potential U{r). For small values of Py, we note that the original step-like 
shape of the potential changes dramatically around p^ = 0. There, it becomes smoother 



0.35 




FIG. 1. Pseudo-potential: a) External potential U{r); b) [W(r, p)]_|__|_; c) [Z//'(r, p)]_| ; d) 

[lY'(r,p)]__. }leTepy/h= W'^ nrn'^ . 



and enlarges the spatial region where the gradient of the pseudo-potential (representing an 
effective electric field) differs from zero. This can be seen as a strong non-locality of the 
potential (or equivalently of the electric field) that is a peculiarity of the graphene band 
structure and reflects the property that a particle around the Dirac point is quite immune 
to the localization effects. In our formalism, we describe this behavior by the presence of 
an effective potential W that becomes more and more non-local when |p| goes to zero. This 
explains way, differing from the scattering process in normal metals, the transport of Dirac 
fermions in graphene sheets is only weakly affected by external sources of disorder (defects 
or impurities). The behavior of the pseudo-potential W around |p| = can be investigated 
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FIG. 2. \W{r, p)]++ component of the pseudo-potential for different values of the momentum py-. a) 
Py/h = 10~^ nm"-*^; b) Py/h = 10~^ nm~^ c) Py/h = 10~^ nm~^. d) Pseudo-potential [^'(r, p)]++ 
for px =0. 



also analytically. By using that 

l[ 1 + e'^'^o.--K+) e-'"^-+(e-(*'^"--'^-+^ - 1) 
e(a+)0t(a-) = - ^ ^ 

2 \ Q''-<t>a+ (^i — e'-('t'a~-'l'a+)^ 1 e-«(0a- 

where we applied the polar notation a"^ = pa±e^'^^^ = ± ^fx^ + ipy, it is easy to see that 

, f -f/(r) 
lim = ^ ^ 

IpKo \^ ^(r) 

so that the in-band component of the pseudo potential vanishes. We remark that this con- 
sideration should not suggest that in the single-band limit the effect of the potential around 
the Dirac point vanishes and the particles move freely. Observing the equation of motion 
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(14) reveals that the pseudo-potential W is non-local both in position and momentum, so 
that the particle motion is conditioned by the values of W in an interval of the momentum 
axes and not just at a point. In particular, only for spatially uniform Wigner distribu- 
tion functions the pseudo-potential becomes local-in-momentum and its single band effect 
vanishes. 

The principal aim of this contribution is to study the effect of the band-to-band transition 
on the stationary current in a graphene sheet. The full quantum mechanical description of 
motion consist of a rather complex set of coupled equations, where a simple interpretation 
of the dynamics is hampered by the presence of the highly non-local operators. In order to 
get more physical insight and to profit of the close analogy between the classical mechanics 
and the Wigner formalism, we consider the so-called gradient expansion procedure. We thus 



expand the functions A and W in Eqs. (13)-(14) with respect to h and limit ourselves to the 



leading order. The study of the full quantum Wigner transport will be addressed in sec. IVl 



If we expand the function A up to the first order in /i, Eq. (13) simplifies to 



[A(p),5l = [A, 5'] - I {VpA, Vr5'} + o{h% (17) 

where curly brackets denote the anti-commutator. In the hypothesis that the external elec- 
tric potential U{v) is regular, we have 

[W, S\ = ihV^U ■ Vp5' + ^ [[9, VpG ■ Vrf/] , S'] + o{h^) . (18) 

This approximation is justified in the limit where the external electric potential U{y) can be 
considered as a sufficiently smooth function, so that only the first-order terms (proportional 



to the electric field) play a significant role in the dynamics. In Eq. (18) it is easy to identify 
the first term with the usual force operator. The second term takes the main quantum 
correction to the classical equation of motion into account. In the following, we will describe 
its physical meaning in terms of the band-to-band transition and exploit its connection 



with the adiabatic Berry phase approximation. The equations of motion (10) become (the 



components of S' are defined in Eq. (16)) 



^ ■ Vr/^ + Vrf/ ■ Vp/± ± z (Sf - S/A , (19) 

dt v^l + ^-2 IpI V / 

dp — 

^ = ^A^ + Vrf/ ■ Vpf + iB (/+ - /-) , (20) 
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where overbar means conjugation and 

e=^^(pAV,C/). . (22) 

^ IPI 

Here, OJt^(^) = /l ± -y^=, ^ = ^^4^ and (p A VrU) denotes the out-of-plane component 
(2;-coordinate) of the vector (p A VrU). 



Equations (19)-(20) extend the semi-classical evolution of a two-particle system in a 
quantum mechanical context. In particular, in the limit of a vanishing electric field, the 
equations decouple and, as expected from a physical point of view, the particle system 
is described by two semi-classical equations of motion. This simple limit eases to attach a 
physical meaning to the various components of the solution. In particular, (/~) represents 
the electron (hole) distribution function in the presence of an external electric field that 
modifies the crystal periodic potential (more precisely, they are the components of the 
Wigner function in a basis, where the two-band Hamiltonian is locally diagonal in the 
momentum and in the position space up to the first order in h). 

To appreciate the advantage of using our quasi-diagonal formalism, we compare the equa- 



tions of motion (19)-(20) with the evolution equations obtained by a direct application of 
the standard two-band Wigner formalism. The Wigner function for a multiband system is 
usually defined as 

where \1/ = (^^) is the two component Schrodinger wave function satisfying = T-L^ . 

Equation ( [23| ) is a straightforward extension of the single band Wigner function, where the 
Wigner transformation is applied componentwise to the density matrix. Up to the first order 
in the equation of motion for the two-component Wigner function writes 

= -f Vr/O + (Vrf/ • V,) + A p , (24) 

'^ = Vrt/-Vp/o-f divf^ (25) 

where we defined the vector = (23? {fai} , 2S {fai} , fn - ^2), fo = fii + f22 and S (3?) 
denotes the imaginary (real) part. The formulation of the two-band Wigner approach given 



in Eqs. (24)- (25) is characterized by the presence of high oscillating regimes. The direct 
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numerical treatment of Eqs. (24)-(25) reveals itself to be a very difficult task. The quantum 
mechanical two-band motion is essentially a two-scale process characterized by band-to-band 
transitions (whose frequency is proportional to the energy difference between states localized 
in the upper and the lower Dirac cones) and the intraband motion of the electrons (that, 
with respect to the tunneling processes, can be considered as a slow dynamical process). 
Furthermore, in this formulation the analogy with the semi-classical evolution of the sys- 
tem (characterized by two uncoupled Liouville equations, one for the particle distribution 
function in the upper cone, and one for the hole distribution function in the lower cone) 
is completely lost. Here, a description of the dynamics where we can associate a certain 
quasi-distribution function to the particle and a different quasi-distribution function to the 
holes, does not apply. One of the most remarkable advantages of the single-band Wigner 
formulation of the quantum mechanics (and was the main reason for which this formulation 
has been introduced) is that in this framework the classical limit /i — t- is easily evaluated. 



As shown by the Eqs. (24)- (25), this is no longer true in the many-band case, where the 
limit h ^ is completely non-trivial. This is due to the presence of the last term of Eq. 



(24). When h goes to zero, the various components of / become more and more coupled 
and the system becomes ill defined. This simple consideration suggests to use instead of 
(f'^, /°), some new unknowns behaving regularly in the limit /i — )■ 0. This can be obtained 
by the partial diagonalization procedure described in this section. 



A. Berry connection in the quantum phase space 

In a crystal where the effective Hamiltonian is expressed by a partially diagonalized 
basis (such as for example graphene or Kane-Luttinger kp models for semiconductors), the 
velocity operator has off-diagonal elements and the electric field mixes the bands, so that 
the expectation value of the velocity acquires an additional term proportional to the field 
and the usual definition of group velocity does no longer apply. The theory of Berry phases 
offers an elegant explanation of this effect in terms of the intrinsic curvature of the perturbed 
band. Furthermore, the Berry connection plays an important role in spin dynamics and in 
describing spin-orbit interactions. We discuss how it is possible to characterize the Berry 
phase in graphene (which is usually studied at the Schrodinger level) by using our kinetic 
description of the quantum dynamics. The formal analogy between spin and band degree of 
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freedom suggests that we investigate the effects of including the Berry phase in the evolution 
of a many-band electron system. We apply the Berry approach to the Hamiltonian symbol 
"H (which is a simple matrix where p plays the role of the adiabatic variable). Berry's 
adiabatic theory states that, if a system is initially described by a certain eigenvector Ui{p) 
of 'H(p), the vector state of the system at time t is given by 

^(t) = Ui{p{t)) e^^^W-s/o -ap(t')] dt' ^ (26) 

where the term 7j is named dynamical phase factor and can be obtained as the path in- 
tegral along the p-trajectory, 7i = / Ajj(p) ■ dp, of the Berry connection A(p) given by 
Ajj(p) = i (Mj(p)|VpUj(p)). In our case, by construction, the distribution functions and 
/~, respectively, are the Wigner functions related to the p-dependent Floquet projectors 



|m+(p)) ( M+(p)| and |u_(p)) ( M_(p)|. Since from Eq. (26) we have that 



\m){m\ = \^^^m]){uMt)]\ 

for these functions, the Berry phases cancel out. On the contrary, the function /* is related to 
the "band transition" operator [^^(p)) ( M-(p)| that, for a given trajectory p(t), cumulates 
a Berry phase equal to 

^ (A.. - A__) - £±M^ - fp A ^1 - ^ IPI ^ A, (27) 



where ^ = VrU. We see that the Berry phase coincides with the "natural" oscillation 
frequency of /*. Our method is thus particularly suited to highlight the role of the Berry 
phase in the evolution of the system. A well known characteristics of the Berry connection 
is the divergence in the proximity of points where the bands intersect. In gapless graphene, 
such a divergence can be found in the neighborhood of the Dirac point p = 0. For that 



reason, from Eq. (21) we see that the natural oscillation frequency A of /* behaves like 
l/|p| when A = 0. 



III. SIMULATION OF GRAPHENE (UNIFORM ELECTRIC FIELD) 

The most common configuration to perform experiments with graphene is constituted by 
the single-layer graphene field effect transistor (FET) Graphene FETs are fabricated 
by standard lithography and a degenerately doped silicon substrate is used to tune the 2D 
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carrier density in the proximity of the Dirac point. We apply our model in the approximation 



of a quasi-uniform electric field (constituted by Eqs. ( 19 )-(20 )) in order to study the quantum 
corrections to the ballistic charge motion in an intrinsic graphene sheet (for which A = 



in Eqs. (19)-(20)) induced by an applied external potential. We prescribe boundary 
conditions in correspondence to the metallic contacts. The contacts are considered as perfect 
charge reservoirs, where the number of particles entering the device are given by the thermal 
equilibrium distribution. By identifying /"*" and /~ with the electron distribution functions 
in the upper (S^) and lower (S^) cone, respectively, we fix their incoming values at the 
boundaries of the simulation domain equal to the Fermi distribution function. Vanishing 
boundary conditions are assigned to the interband function /*. 



We consider a simple device consisting of a graphene sheet suspended by two ohmic 
contacts at a distance of 1 fim. The bias voltage U is applied between the contacts. This 
prototype of devices has been experimentally probed in [9J. The presence of interfacial 
phonons in the substrate reveals itself to be an important source of limitation for the charge 
mobility in graphene. However, suspended graphene offers the considerable advantage that 
the interactions between the underlying substrate and the graphene sheet are completely 
eliminated. Up-to-date lithographic technique allows the fabrication of high quality graphene 
sheets suspended on a silicon substrate where the mean distance between the flat graphene 
sheet and the substrate is around 150 nm. Under this conditions, we can safely assume that 
no phonons are transmitted to the graphene sheet from the substrate. At room temperature, 
mobilities of suspended graphene are close to 10"^ cm^V^^s^^ and are limited by acoustic 
phonon scattering. Mobilities of such an order of magnitude imply that electrons can travel 
from one contact to the other by suffering only a few scattering events. This evidence 
justifies the study of ballistic transport in suspended graphene. As a further simplification, 
we assume that the particles move under the action of an external electric field S directed 
along the x direction and independent of the y variable. In this case, we can assume §^ = 
which greatly reduces the numerical complexity of the system. The equations of motion 
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FIG. 3. Stationary solution for graphene for an applied potential Vq = 0.3 eV. Snapshots of the 
(upper plot), 1 — /~ (lower plot) distributions on the x — px plane, for Py/h = 0.1 nm^^. In the 
left plot we represent the contour lines and in the right plot the 3D representation of the solutions. 



(19)-(20) simplify to 



9t ^ ^vi + vi pI + pI \ ^pI + p'l 



dt dpx 2(p2+p2)3/2^ 



(28) 



(29) 



where 



and p = (pxiPy), I" = {x,y)- In the figs. [3]|4]we depict the stationary values of the electron 
distribution and the hole distribution (1 — /^) for an external applied potential Vq = 0.3 



Px ~^ Py "' 
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FIG. 4. Stationary solution for graphene under the external potential Vq = 0.3 eV. Snapshots of 
the (a-b) and 1 — (c-d) distributions on the x—p^ plane for Py/h = 2- 10^^ nm^^. In the left 
plots we represent the contour lines and in the right plots the 3D representation of the solutions. 
In (e-f) we plot the semiclassical solution (without band transition). 



eV for different values of the parallel momentum py. In the left plots we represent the contour 
lines and in the right plots the 3D representation of the solutions. In our simulation, the 
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lattice temperature T = 300 K. It can be clearly seen that the electrons split in two parts: 
the particles with velocities parallel to the electric field which are accelerated, and anti- 
parallel ones which are refiected back by the potential barrier. Further, due to the presence 
of interband Klein tunneling, also interband particle transitions between the bands and 
S~ are possible. Since the relation between the velocity and the momentum for a hole is 
the inverse of that for an electron during this interband transition, the momentum parallel 
to the barrier is conserved and the velocity of the quasiparticle is inverted. Due to the 
larger number of particles in the lower cone, we observe a net fiux of particles from E~ 
towards As expected, interband transitions become a dominant phenomenon around 
the Dirac point p = 0. In fact, in correspondence to high values of Py (depicted in figure 
|3|, the distribution functions look very similar to their classical counterparts and quantum 
corrections are negligible. On the contrary, for smaller values of Py (see fig. |4| a fiux of 
particles from the S~ band to the band is clearly visible. To highlight the effect of 
Klein tunneling, in fig. |4]-e-f we present the distribution functions of electrons and holes 
under the same condition as in fig. |4]-a-b but in the semi-classical approximation (without 
tunneling). One of the advantages of our approach is that now the Klein tunneling effect 
can be described by the familiar language of classical mechanics. In fact, from fig. |4]-c-d 
we see that, in order to overcome the potential barrier applied between the two contacts, 
a large number of particles belonging to the S~ band leave this band and a corresponding 
increase of the related hole distribution function (1 — /~) is observed. These particles are now 
accelerated by the same electric field in the final part of the device {x = L) and contribute 
to increase the particle distribution /+. 

A. Current and density 

The density and current of particles in the upper (lower) band, denoted by n+ (n^) and 
j+ (j^), respectively, can be obtained from the Wigner functions as 




(31) 



(32) 
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FIG. 5. a-c): Particle density in the S"*" band (n^) for different applied potentials, d): 3D plot of 
n+ as a function of the applied potential and the position. 



The continuity equation for the charge can be deduced from the system of Eqs. (19)-(20): 

V.fT A^t/idp, (33) 



dt 



-M[f] = ^(pAVrf/),53{/V^p} . (34) 

In particular, corresponding to a stationary solution (where = 0), the total current 
j* = j+ + becomes uniform (Vrj* = 0). Figure [s] shows the stationary charge density 
profile n"'"(x) in the intrinsic graphene for different applied voltages. For low voltages, the 
behavior of the particle density is essentially semi-classical: with the increase of the external 
field, the electrons cumulate near the source contact and a charge depletion in the channel is 
observed. The nearly total depletion of the drain contact is reached for an applied potential 
of 0.15 eV (snapshot a of fig. [s]). In correspondence to a further increase of the applied 
potential, the quantum Klein effect starts to play a relevant role in the shape settlement of 

20 



the stationary density profile. In particular, for a potential greater than 0.25 eV (snapshot 
c of figure [s]), we observe a monotone increase of the charge density inside the channel. 
This effect is due to the particles, initially localized in the S~ band, that are injected in 
the upper cone as a response to such a strong electric field. We see that in this regime 
of a strong external potential, the quantum correction to the density becomes comparable 
with the total charge present in the device. Finally, in fig. |5]-6 we highlight the presence of 
density oscillations in the proximity of the ohmic contact. 




We focus now our attention on the quantum correction to the total stationary current. 
It is well know that, when scattering processes are neglected, no steady state can exist in 
the graphene bulk. If a uniform electric field is applied to a spatially infinite sheet, the 
momentum of the particles would increase indefinitely and the current would show Bloch 
oscillations [21] . Nevertheless, when the real band structure of the graphene is approximated 
by an unbounded bi-conical shape, the saturation of ballistic current is reached. In fact, no 
matter how much they are accelerated, particles produce always the same amount of current. 
The upper limit of the current is obtained when all the particles entering the device through 
the source contact reach the drain and, at the same time, the drain incoming particles are 
reflected by the barrier. The current {Jgat) and density {riin) related to the incoming particles 

21 




0.1 0.2 0.3 0.4 

Potential (eV) 



FIG. 7. I-V Characteristic: Comparison of the quantum solution (continuous green hne) with the 
semi-classical solution (dashed blue line). 



distribution /j„ at the source contact are 



-1 



^» = ^^ J |p|/.n(|p|,^)d^d|p| = ^(^)Yp[l + e''^'"""'^-"^)j dp, (35) 

•J sat — ^in- {oO) 
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In figure |6] we plot the saturation current Jgat versus the chemical potential fi for a tem- 
perature of 300 K. The ballistic saturation current shows a considerable increasing when 
the chemical potential is augmented. On the contrary, numerical simulations proof that the 
quantum correction to the total current induced by interband tunneling is almost insensitive 
to a variation of the chemical potential. This can be understood if we note that the Klein 
tunneling in presence of a (almost uniform) slowly varying electric field concerns particles 
whose energy is located around the Dirac point (or equivalently, particles whose momentum 
is nearly zero). For low temperature, if the chemical potential is above the Dirac point, the 
number of such particles is almost independent from fi. In this contribution, we will focus 
our attention to a quasi-intrinsic graphene sheet, for which the interband current is of the 
same order of magnitude as the saturation current. 

The stationary I — V characteristic of the device (intrinsic graphene) at the temperature 
of 300 K is depicted in fig. [7], where the current / fiowing through the device is plotted as 
a function of the bias voltage Vq applied between the source and the drain contacts. We 
compare the solution of our quantum system with the classical motion. For this purpose, 
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we plot the results obtained by discarding the interband transitions (dashed blue line) and 
including the multiband corrections (continuous green line). Our simulations show the 
importance of including the tunneling process in a realistic simulation of the current in 
intrinsic graphene. In particular, in the case of quantum transport, no saturation is observed 
and the current grows with increasing external potential. 

B. Non intrinsic graphene 




FIG. 8. Stationary solution for an applied potential of Vq = 0.3 eV and ^ = 0. Snapshot of the 
/"•" distributions in the plane px — Py, at different positions x along the device: (a) x = L (source 
contact) (b) x = L/3 (c) x = —L/3 (d) x = —L (drain contact). 



In order to give a clearer description of the two-band motion and to compare the solution 
of intrinsic graphene with doped graphene, we represent the solutions in the Px—Py plane. 
In particular, in fig. |8]we represent the stationary solution for the distribution function /"*" 
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FIG. 9. Stationary solution for graphene under the external potential Vq = 0.3 eV and /i = 0.6 eV. 
The snapshots represent for different positions x along the device : (a) x = L (source contact) 
(b) X = L/3 (c) X = —L/3 (d) x = —L (drain contact). 

in intrinsic graphene under the same condition as in the figs. |3]|4] (the external potential 
Vq = 0.3 eV). The snapshots show the /"*" distribution at different positions x along the 
device: (a) x = L (source contact), (b) x = L/3, (c) x = —L/3, (d) x = —L (drain contact). 
We see that the particles entering the device from the source contact are accelerated by 
the potential and leave the device at x = —L without reflection. On the contrary, particles 
injected in the graphene sheet from the drain contact have not enough energy to overcome 
the potential barrier and are reflected. Based on these general considerations, we see that in 
the snapshot of fig. |8]-a) only one electron beam is visible. The following cuts of the solution 
along the x-plane, toward the drain contact (fig. [8]-b,c), show that a new particle beam 
around p = appears. This second pulse describes the particles coming from the source 
contact and cumulate along the channel. Similar consideration hold for the S~ particles (we 
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do not report here the E~ distribution /~), with the difference that in this case the band is 
almost full and, for small values of the momentum, a large number of particles are able to 
overcome the small potential gap between the E~ and the S"*" band (which is equal to 2f^|p|). 
As a consequence, S~ particles coming from the source contact, x = L, are decelerated by 
the electric field and instead to be completely reflected back to the source contact, they 
leave the S~ band. These particles, now belonging to the S"*" band, are accelerated by the 
electric field and contribute to the increase of the high energy electron beam depicted in fig. 
[sj-d. The opposite Klein processes, where the particle flow is directed from the S+ to the 
S~ band is still present but with a smaller intensity and not visible in our plot scale. In 
fig. [9] we represent the stationary solution for the distribution function in non-intrinsic 
graphene for an external potential Vq = 0.3 eV and for fi = 0.6 eV. As shown in fig. |6j 
for higher values of the chemical potential /i, the semi-classical intraband current quickly 
increases. In contrast to the intrinsic case, now we observe that particles in the S"*" band 
populate higher momentum levels and quantum tunneling becomes less significant. 



IV. EFFECTIVE MODEL 



A direct solution of the system of Eqs. (19)-(20) and its application to electron transport 



in a graphene sheet demands a high computational effort. In this section, by investigating the 
general properties of the solution, we derive some asymptotic limits where an approximated 
version of the equation of motion applies. The major problem arises from the approximation 
of the equation of motion for the interband function The diagonal functions (which 
can be considered as a straightforward generalization of the distribution function of electrons 
in the upper and in the lower part of the Dirac cone) share similar properties with their 
classical counterparts, and are rather smooth and stable. On the contrary, the interband 
function /* shows high frequency oscillation regimes. 

In the study of the electric properties of a solid, like the current-voltage I — V charac- 
teristic and the conductivity, it is often of primary interest to obtain a correct description 
of the non-equilibrium stationary state reached by the system in response to an external 
perturbation field. In the case of the I — V characteristic, the external perturbation is rep- 
resented by the gradient of the applied potential. The knowledge of the stationary I — V 
characteristic is crucial for engineering applications of a material and to its integration into 
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a network. In particular, different approaches should be adopted if a system is characterized 
by a single time scale according to which all the interesting observables evolve, or if some 
observables evolve much faster than the others. In the latter case, these variables identify 
some "internal dynamics" of a multi-scale process. In our system, the interband function /* 
is a strongly oscillating function and its "natural" oscillation frequency A depends on the 
momentum p. This reflects the general principle of quantum mechanics that a wave function 
containing a superposition of states with different energies, oscillates with a frequency which 
is proportional to this internal energy difference. In our case, /*(r, p) describes a mixture of 
states belonging to the upper and the lower cone. At a given position r, their mean energy 



difference is equal to 2i;i?|p|. This appears explicitly in Eq. (21). Because of the high value 
of the Fermi velocity in graphene, this term induces a dynamical evolution of that can 
be considered to be considerably faster than the other processes induced by the external 
field (we remark that the identification of the different time scales in which the two-band 
quantum system evolves, is practically infeasible with the usual definition of many-band 



Wigner functions given in Eq. (23)). Since describe states with similar energy, in view 
of Eq. (28), no "natural" oscillation frequency is present in the equation for 



We are interested in deriving an approximated formula that integrates the function 



/j(r, p,t). Equation (29) can be recast in integral form as 



f(r,p,t)= r '\^^o^^P^+£r,Py)dr^(^p^^st')f{r;p, + St\py;t-t')dt', (37) 

^0 



where A is given by Eq. (30), = — f and 



After some algebra we obtain 



(38) 



i _ 

K ' ^ Pa ' ' Py 

obtain an asymptotic expression for the function /* in the limit t — t- oo. A simple analysis 



with a = /3 = £31, 7 = and ^{(3) = (3^/lT]3^ + ln(/3 + y^lT^). We intend to 



of Eq. (38) reveals that the function /* displays two qualitatively different behaviors if Py 
is grater or smaller as a certain value A. In particular, /* is smooth if < A and becomes 
strongly oscillatory otherwise. In the following, we make this statement more precise. In 
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FIG. 10. Schematic representation of the region 



the hypothesis of smooth the long-time behavior of the function /* can be estimated by 
studying the integral 



This expression suggest to estimate X by means of the stationary phase approximation. This 
approximation applies when ^ 1 and the exponential is fast oscillating in the scale 

of the polynomial decay {u^ + Explicitly, this condition gives \py\ ^ A = or 

\px\ ^ A'^/py. In order to give an analytical estimation of we divide the plane p^ —Py in 
the interior and exterior part of the region VL defined by 

r]={p:b,|>A or > AVb.l}, (40) 



which is depicted in fig. 10 For p G the stationary phase approximation applies and 



the integral of Eq. (38) can be easily estimated. On the contrary, for p ^ a different 



approximation is adopted. Some numerical tests are presented in Appendix |VIIB where the 



validity of the approximation procedure used for the derivation of the asymptotic evolution 
equation is investigated. 

The momentum p^ is evaluated along the trajectory Pxif) = Pxito) +£{t — to). For p G f2 
the phase velocity A changes in time and increases indefinitely for t going to infinity. The 



form of Eq. (38) reveals that the long time-behavior of the solution /* is dominated by 



the exponential term. By expanding the function in the exponential up to the third order 
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around the stationary point m = 0, we obtain the following approximation for 

Py + ^Px 3/2 ^. (a/7;^\ „ +^ j%m 



r(r,p,t) ^ 7re(-fp^)sgn(p^pj^) 



a 



Ai (^v^) /'^ (r,p,t) e't' 



\/pI + 

where ^ denotes the Heaviside step function and Ai(x) the Airy function: 



(41) 



27r Ai ix] 



dt 



(42) 



Up to the first order of the electric field, / in Eq. (38) can be approximated by /q 



fo — fo , where Z,^ = (l + e^^^'^') are the Fermi distributions. By expanding /q around 
the stationary point p^: = 0, we obtain 



foir;Px + £t,Py;t) 



Px=0 



4 \Py 



(43) 



For a temperature T = 300 K, an electric field £^ = 0.1 eV/im"^ and a parallel momentum 
Py/h = 0.1 nm~^ (which are the typical values for graphene), the previous equation reveals 
that, around the stationary point px = 0, the /q function evolves in a time scale of picosec- 
onds. This time scale is considerably smaller than the "natural" frequency A, which is of the 
order of femtoseconds. These considerations suggest to simplify the evolution of the system 
by assuming that the time evolution of the diagonal functions is smooth compared to 
the time evolution of In this hypothesis, we consider an asymptotic model where the 
function is assumed to be constant around the stationary point Px = 0. 

For p ^ Q (for the sake of simplicity we assume also |p| — t- 0), we approximate the 



function A given in Eq. (30) with the dominant contribution 

Py 



A--S- 



pI + pI 



(44) 



and obtain from Eq. (38) after some algebra 



/*(r,P,i) 



TT 



sgYi{£py) — tan ^ 



Px 

Py 



2 ^pI + pI 

This approximations, together with Eq. (41), lead to the following equations of motion 



(45) 



dt 



T{Px,Py) 



Px 



^pI + pI 9x dpx 



r 



£ Py 



|sgn(£py) — tan^ 



^1 I Px 



-^^{fi^^ (^) 0{-Spx)sgn{pxPyW {e^^^(«} pen. 



(46) 
(47) 



28 



In contrast to the full quantum mechanical formulation of the dynamics given by Eqs. (19)- 



(20) or Eqs. (24)-(25), this form of the approximated equations of motion reveals itself 
to be a simpler and easily understandable description of the interband coherent quantum 
tunneling phenomena. Here, the transition of a particle between the two bands is modeled 
by a balance equation, where the transition probability is given in terms of the "tunneling 
scattering rate" T. The scattering processes is described in a simple way; in the presence 
of an electric field £ (directed for simplicity along the direction x), the component of the 
momentum parallel to £ changes according to the Newton law Pxif) = Px(^o)+^(^~^o) and 



when — tunneling occurs. Basically, in Eq. (47), we distinguish between small {py < A) 
and large {py > A) parallel momenta (we remark that py is unaffected by the presence of 
the electric field and plays the role of a parameter). As explained below in more details, for 
Py < A the transition rate T is proportional to l/|p|. Under this condition the band-to-band 
transition becomes highly favorable and can be interpreted by a quasi-instantaneous process 
taking place when the particle is at rest. In the opposite limit Py > A, a. complex pattern of 
interference between the two bands appears, which gives rise to the highly oscillatory shape 
of This oscillatory behavior is captured by the stationary phase approximation, and 



appears in Eq. (41) through the phase a^{l3) which modulates the transition rate T. This 
part of the solution is the origin of some numerical noise that can be observed in a direct 



numerical discretization of the equation of motion (19)- (20). In particular, concerning the 
application of these equations to a graphene sheet by using a reasonable size of the mesh grid 
(the number of the grid points for the p axis being of the order of 10^ — 10'^), the function 
Q;^(/3)/(27r) covers many periods within each p^ discretization cell, making the numerical 
solution quite inaccurate. This analysis suggests that a possible solution to this problem is 
to substitute the function /* by its Gaussian convolution around each discretization point 
in the Px — Py plane. Anyway, despite the uneasy form of the strong oscillation regime 
prevents any interesting phenomena to emerge at the macroscopic scale (the expectation 
values of any observable being expressed by the integral of the Wigner functions, so that high 
oscillating contributions average to zero). Moreover, for py ^ A the transition rate T, and 
consequently the interband tunneling probability, decreases exponentially. This exponential 
decay agrees with the well known Landau- Zener formula for which the transition probability 



Vppy 



is proportional to e if [22]. To go more into details, the formula of Eq. (47) approaches 
the Landau-Zener probability in the limit of \py\ going to infinity. In this case, T can be 



29 



simplified by using the asymptotic expression Ai(x) ~ %^^;;t7I"- We obtain 

c I — \ 3 he 

nP.,Py) - -^IT^V r 23/4 ^ g(-gp.)sgnto)g?{e-t^(^)} , 

where we note that the exponential decay is well represented, but with a slightly different 
rate (we found instead vr). 

The interesting regime for studying the Klein tunneling process is given by |p| < A. From 



Eq. (47) we observe the emergency of some interesting limits revealing deeper insight into 



the physical description of the tunneling processes provided by our formalism. We consider 



Eq. (47) for p.j. = 0. In this case, the transition interband probability T becomes 

which goes to infinity when py approaches zero. We show that this divergence reflects the 
well know property that a particle, whose trajectory passes exactly through the point in 
the energy spectrum where the upper and the lower cones touch (F point), has a unitary 
probability to pass from one band to the other one. At the F point, the distinction between 
the upper and the lower band becomes artificial. For this reason, the distribution function 
for S+ particles and /" for particles should be equal at p = (they represent the 
same quantity). Any configuration of the system where f~^{r, 0) ^ /^(r, 0) is unphysical. In 
our formalism, it is easy to see that the equation of motion ensures automatically that this 
condition is satisfied at any time. In this context, we can think about our two-band model 
as a system representing the evolution equation for two populations of strongly interacting 
particles. The scattering kernel is now written in the relaxation time approximation, where 
the relaxation time goes to zero. The equation of motion for the difference — is 

In the limit T — t- oo (Chap man- Enskog limit), we neglect the drift term and the previous 
equation gives 

and up to the order o(l/T) we obtain /"*" = A careful analysis of the origin of the 
divergence in the transition rate T reveals that the divergent term is exactly the Barry 

30 



connection associated with our two-band system (see section II A). This can be seen by noting 



that in the approximation of the phase A given in Eq. (44), we retain the Berry connection 



1^ (p A Vrf/)^- We discard the contribution ^|p| that equals the difference of energy 
between states with same momentum but localized in different bands. We remark that 
even in our model a divergent term appears explicitly. This does not lead to an unphysical 
result that is usually found in similar situations. Our model is able to include explicitly the 
divergence of the Berry connection by simply forcing the solution to be equal to zero in 
the region where the Berry phase is not defined. 

We consider now more generally the limit Py — ?■ 0. Care have to be taken in evaluating 
this limit. Instead of considering directly the function T, it is convenient to consider the 
main integral value of this function in an interval = [px — e, Px — e] and we let e go to zero 
at the end of the calculation: 



T{Px,Py) 



2e 



Px+t 



T{Px,Py) dpx 



1 £ 



2e2 



-L 2 

sgn{£py) ta.n^^{v) — - (tan^"^) (t>) 



(Px-e)/py 



If Px 7^ we choose e < \px\ and obtain 



lim T{px,Py) = 0. 

Py^O 



For Px = the limit yields 



e lim T{px,Py) 

Py^O 



1^1. 



These considerations show that the correct limit for the transition probability T is given by 



\im T{px,Py) 

Py^O 



TT 



imp.) ■ 



This form of the transition rate put in evidence that the band transition around the Dirac 
point is a strongly localized process. Based on the previous analysis of the behavior of the 
transition rate T, we further simplify its expression by evaluating the main value of T. We 
integrate the transition probability with respect to Py, 



T{px,Py) dpy 





s 


TT 


\Px 


2 



[1 - sgn{Spx)] log 



A 

Px 



where, for the sake of simplicity, we approximate tan ^ ( ^ 



^ ^sgn{pxPy). We note that 
T is non zero only if the momentum of the particle Px has the opposite sign with respect to 
the electric field S. According to macroscopical considerations, this means that a particle 
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undergoes a transition only if it is decelerated by the field (the transition takes place when 
the momentum of the particle can be considered to be small compared to A). We remark 
that this consideration applies irrespective of the cone in which the particle belongs to, since 
the classical equation of motion for the momentum is p = i^(x) in both cases. 



V. FULL QUANTUM SOLUTION 



We focus our attention to the numerical solution of the full quantum mechanical electron- 



hole pair evolution. We consider Eq. ( 10 ) without any further approximation. For the sake 



of clarity, we report here the equations of motion 

dS' 



(48) 



where 



n' =U'{r,p) + A{p) , (49) 
W (r> P) = 7^ / (P + et (p - ^f,^ U{v')e^^-^'y^ dM dr' , (50) 



and 



1 



(2 



( r - ^77, p + j S' (r, p ) - S' (r , p') H' (t + ^r?, P - 



X 



gi(r-r')-M+i(p-p') r? _ 



In order to obtain a numerically tractable model, instead to solve the full 4-dimensional 
system (two-dimension both in position and in momentum), we consider a simpler case where 
the solution is uniform along the y direction (but non-constant with respect the momentum 
along the same direction) and we solve the reduced system in M? . In this hypothesis, we 
have 

S' (r,p) = S' ir^;p^,py) , 
W (r,p) = U' {r^;p^,py) 

and the equation of motion simplifies to 
1 



W [r- ^ri,p+^fx,pyj S'{ii,ri,py) - S'{ii,ri,Py)U' (r+ ^rj^p- ^fx,Py 



xe'^'+'P" d/i dr] 



S'if^,v,Py) = / S'{Cp'^,Py)e-'^^>^^'^^^ dr^ dp',. 
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with 

U\r,p,py) = ^ I e + ^ii,p)j et (^p - ^^,,p)j j U{ry^^-^> dr' dfi , 

where we put in evidence that in this case the coordinate Py plays the role of a parameter. 
An analysis of these expressions reveals that in order to obtain an efficient numerical scheme, 
it is convenient to impose 

Ap = pA^/2 , 

where a, /3 are integers (or inverse of integer) and denotes the size of the numerical 
discretization of the z axes. Since the x — fi and p — t] are conjugate variables, the discrete 
Fourier transform (DFT) requires (see for example 
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Ar 



NrA^ ' 



A 



and we obtain that the following relationship should be fulfilled: 

Np a ' 

Finally, we consider the following first-order (in time) solution of Eq. (48) that shows itself 
to be particular stable and weakly affected by numerical noise: 

S' (r,p;t + A,) = /e-^"'('-t'^'^+t'^)5'(/x,r/,t)e^"'(^+t^'P-|M)e-/^+^f. d/i d.^, 

(27r) J 

where the matrix e^^'^'^'^^ is evaluated by the formula 

^'^^'(''P^ - I cos(|u|)Mo + ^sin(|u|)^^ e*"° , 



|U| 

u = Y tr{W'(r,p) - a} , 
«o = ^ tr{W' {r,p)} 

with tr denoting the trace of the 2x2 matrix. 

In order to present the structure and to give a general impression of the form of the 
full quantum mechanical solution based on the matrix Wigner function S', we consider 
a standard text-book case, in which a minimum uncertainty Gaussian packet impacts a 
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FIG. 11. Contour plot of /+ (left plot) and /" 
t = 700 fs; d) t = 900 fs. Here py/h = lO^^ 



' (right plot) for a) t = 300 fs; b) t = 500 fs; c) 
1 and Vq = 0.1 eV. 



potential barrier. The study is depicted in fig. 11 In particular, we consider as initial 
condition a Gaussian pulse in the upper graphene cone (S"*" band) localized around the 
position Xo = —40 nm and momentum po/h = 10~^ nm~^ and with a parallel momentum 
Py/h = 10~^ nm~^. Furthermore, we assume a vanishing initial condition for the S~ band 
and a vanishing band-to-band correlation (represented by the function /*). The shape of 
potential barrier used in the simulation is depicted in fig. [T| but here we consider a lower 
barrier of 0.1 eV. The wave is initially localized in the zero potential region {x < 0) and 
is directed against the potential step. The height of the barrier is chosen in order to be 
smaller than the mean kinetic energy of the wave packet, so that it could be overcomed by 



the Gaussian packet. In fig. 11 we display the solution for different times (from sub-panel 



a to d). In particular, we show the contour plot of (on the left side) and / (on the 
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right side). The solution shows that the S"*" packet overcomes, as expected, the potential 
barrier but also generates a transmitted particle beam in the S~ band. We note that the /~ 
function is initially generated very close to the potential barrier, where the main momentum 
of the E"*" particles tends to the minimum. In order to give an impression of the relative 



width of the two solutions, in fig. 12 we depict the 3D version of fig. 11 



In the second numerical test (depicted in fig. 13), we consider a more stressed case 
consisting of a higher barrier (difference of potential equal to 0.3 eV) and a lower parallel 
momentum Py/h = 10""^ nm~^. According to the previous discussion, when the momentum 
of the particles approaches the Dirac point p = 0, the coupling between the functions 
and /~ increases considerably. The emergence of a divergence in the coupling terms entails 
that the numerical solution becomes more and more critical. This consideration explains 
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FIG. 13. Contour plot of /+ (left plot) and /" 
t = 700 fs; d) t = 900 fs. Here py/h = lO"'' nm" 



' (right plot) for a) t = 300 fs; b) t = 500 fs; c) 
1 and Vo = 0.03 eV. 



the emergency of the complex interference pattern that is observed in the phase plane where 
several ripples appear both in the upper and in the lower cone distributions. Anyway, the 
main classical features of the solution are preserved (especially in the classical-like region 
|x| > 50 nm). We see that the incoming particles impact the potential barrier and are 
reflected back. Besides, a transmitted pulse in the lower cone is generated. We note that 
both pulses stay "mainly" positive with some residual oscillation induced by the band-to- 
band interference. 

Finally, we address our attention to a more realistic case, where an ohmic contact is 
localized on the left part of the domain at x = —50 nm. We model the contacts in the 
usual way by assuming an incoming thermal equilibrium distribution for the particles 
(to highlight the Klein phenomenon, we artificially impose vanishing boundary condition 
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FIG. 14. Contour plot of /+ (left plot) and /" (right plot) for a) t = 300 fs; b) t = 500 fs; c) 
t = 700 fs; d) t = 900 fs. Here py/h = 10^^ nm.-^ and Vq = 0.03 eV. 



for the f~ distribution). We put evidence into the effect of the full quantum band-to-band 
tunneling (and also to study transient effects) by initially discarding any band-to-band 
effects. For this purpose, we take as initial condition for /"*" the full quantum single band 
thermal equilibrium distribution. For t > we allow particles to pass from one band to 
the other one by solving the complete two-band system. We observe that, after a transient 
regime, a constant flux of particles is generated in the band. This new particle beam can 
be interpreted in the classical language as a flux of particles initially created around Px = 
and subsequently accelerated by the electric field (that is a barrier for the particles in the 
upper cone and an accelerating field for those localized in the lower cone) and propagate 
afterwards freely in the classical (x > 50 nm) zone. In particular, we remark that the 
solution shows the nice property that the transmitted particle beam stays manly positive 
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FIG. 15. 3D plot of /+ (left plot) and /" (right plot) for a) t = 300 fs; b) t = 500 fs; c) t = 700 
fs; d) t = 900 fs. Here py/h = lO'^ nm-i and Vq = 0.3 eV. 



(at least within the numerical precision of our simulation). Here, the hight of the barrier is 



equal to 0.3 eV and Py/h = 10 ^ nm ^. In fig. 15 we depict the 3D profile of the solution 



VI. CONCLUSION 

In this contribution, the ballistic transport of electrons in graphene by including quantum 
effects is investigated in terms of the Wigner formalism. The resulting formulation reveals 
itself to be particularly close to the classical description of the particle motion. Special atten- 
tion is devoted to model the Klein tunneling and to study the correction to the total current 
in intrinsic graphene induced by this phenomenon. Due to the high numerical complexity 
of the resulting system of equations, an approximated closed-form solution is obtained. The 
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simulations show that for an intrinsic graphene in the presence of a strong electric field, our 
model predicts a non- negligible correction to the charge inside the device. Some numerical 
experiments are performed where the evolution of a Gaussian pulse in the presence of a 
potential barrier is investigated. The numerical solutions show that our formalism shares 
some nice properties with the classical solution like smoothness and positivity in the regions 
sufficiently faraway form the potential barrier. 



VII. APPENDIX 



A. Derivation of Eq. ( 10 ) 



In this section, we derive the equation of motion (10). We recall some properties of the 



Weyl operator algebra. Coherently with the notation used in sec. [11} we will denote by 
^ = W [^] the operator associated with the phase-space function ^(r,p). The following 
property holds true 



W 



-1 



AB 



Ai<B. 



(51) 



In the hypothesis that A and B are sufficiently smooth, the Moyal product defined in Eq. 



(11) admits the following /i-expansion: 

ih 



A^B = Y^ 



r^(r,p) 



i3(r,p) 



(52) 



EE (f )" ^(:)^(r'P) ■ V:)"-^ [% . V:)^.(r,p), (53) 

n fc=0 ^ ^ • \ / 



where the arrows indicate on which operator the gradients act. In particular, if both op- 
erators depend only on one variable (r or p), the Moyal product becomes the ordinary 
product 



^(p)^S(p) = ^(p)S(p). 

The Moyal product can be expressed also in integral form: 
1 ■ ' ^ ^(^v^-^v 



(54) 



Ai^B = ^^A{y,tp) e^y"'"^'"^"'') 
(2vr) 



j B (r', p') e^('--r')-M+i(P-P')-r? (ip' 
j A(^r-^r],p+^^?jB (r', p') e^('-'')-'^+^(P-p')-^ dr' dt] dp' , (55) 
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where we used the expansion of Eq. ( 52 ) in the expression 



In same way we obtain 



Bi<A= [ B (r', p') ^ ( r + ^ry, p - ^fi] e^('--r') M+*(P-P') r? (55^ 

(27r) J V 2 2 J 

We evaluate the Moyal symbol T-L' of the transformed Hamiltonian T-L' = OH©^. By using 



Eq. (51 ) we obtain 



n' (r,p) = e*?{^et, 



(57) 



where Gt = W 



-1 



et 



e 



This can be verified easily 



. In particular, 6^ = ^W^^ 
by applying the Weyl operator to the relationship GO^ = X, where X denotes the identity 
operator. We obtain 



x = e^Gt = e et = e 



(58) 



where in the second equality, we used the expansion of Eq. ( 52 ) . The symbol G defined in 



Eq. \\8« does not depend on the spatial variable r. Eq. (57) thus becomes 



n' (r, p) = ^ [Ho(p) + f/(r)] ^ = A(p) + e ^ f/(r) ^ e^ 



(59) 



where we used Eq. (54) and Eq. ([T]). Proceeding as in Eq. (55) we obtain 

W = e(p) ^ f/(r) ^ et(p) = e(p)e-^(^^)f/(r)e^(^^)0t(p) 

= [ e(p)eKMf/(r')e*(^-^')-'^e-K'^^)et(p) d/^ dr' 

(27r) J 

= ^ / (p + (p - t/(r')e^(^-^') '^ d^ dr', 

where we apphed the identity 

f/(r) = / / f/(r')e^('-'') '^ d^ dr' . 
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FIG. 16. Stationary values of the function /* in the presence of a uniform electric field E for 
Pylfi = 0.2nm~-'^. 



B. Numerical study of the asymptotic model 



We present here some numerical tests that validate the approximations used in sec. IV 



and show some characteristic features displayed by the function /* in correspondence to 



the different limits previously discussed. We solve Eqs. (37) with high numerical precision 
and obtain the stationary solution /* in the presence of a uniform electric field S. We 
assume for simplicity = 1. The numerical results show that, as expected, the function 
/* displays high-frequency oscillations along the Px axis. This behavior becomes more and 
more evident when the parallel momentum py goes to zero. In graphene, the band-to- 
band transition probability approaches one for p = 0. For this reason, small values of Py 
characterize the interesting regime, when we study quantum corrections to the interband 



current. Moreover, Eq. (30) shows that for py going to zero, the function /* oscillates with 
a period of Ap^ oc The monotonic increase of the oscillation frequency along the Px 

axis for increasing values of Px (we recall that the electric field is directed along the x axis 
and we are evaluating the integral along the trajectory Px(t) = Pxiio) + £ii — io)) makes the 
direct numerical approximation of /* quite delicate. Besides, the high oscillating behavior 
of /* contrasts the form of the diagonal functions that stays smooth even in the presence 
of significant band-to-band transitions. Thus, it ensures the validity of the stationary phase 



approximation used in sec. IV 
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FIG. 17. a-b) Numerical solution of /* for Py/h = 2-10 '^nm ^. c-d) Approximated expression of 



/* given in Eqs. (41) and (45). 
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In our simulations S = 0.3 ■ 10^^ eV nm~^ and A = = 2 ■ lO^^nm"^. Figure 

shows the typical form of the solution when p & Q {py = 0.2 nm~^). The function is highly 
oscillating with a period of nearly 10~^ nm~^ which is considerably small with respect the 
typical spatial variation of the classical distribution functions. To make a comparison, the 
functions in graphene at the temperature of 300 K vary on a resolution scale of the 



order of 10 ^nm ^ (see for example fig. pi). According to Eq. (41), the numerical solution 



confirms the exponential decay of /* for increasing values of Py. In fig. 17 we represent the 



solution for Py/h = 2-10 ^nm ^. In particular, in fig. 17-a-b we display the function / 



evaluated numerically (figure 17 -a is the snapshot of the zoom of fig. 17 -b in the central 



region) and in fig. 17-c-d we depict the approximation of /* obtained by Eqs. (41)-(45). A 
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glance to fig. 10 shows that if Py < A (as in the present case), the small values of px does 



not belong to Q and should be approximated by Eq. (45). On the contrary, for increasing 



values of Px, Eq. (41) applies. For the sake of clearness, we marked in fig. 17 the boundary 
of the region Q. We remark that around px = 0, where px — Py (and thus p ^ Q) the 



expression of /* is well reproduced by a simple pole (Eq. (45)). For increasing values of 



the function starts to oscillate in same way as described by Eq. (41). 
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